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INTRODUCTION 


Localized buckling of a delaminated group of plies creates high interlaminar 
stresses around the boundary of the delamination. This often causes rapid, extensive 
delamination growth. A measure of the intensity of the interlaminar stresses at the 
delamination boundary is needed to predict delamination growth. Since the stresses 
at the delamination boundary are singular (at least mathematically), calculated 
stresses there have little meaning. Strain-energy release rates are finite param- 
eters which characterize the intensity of the stresses near the singularity. Hence, 
many of the current efforts to predict instability-related delamination growth are 
based on strain-energy release rates. 

The problem of instability-related delamination growth has been studied both ex- 
perimentally [eg. refs* 1-4] and analytically [eg. refs. 3-7]. The predominant con- 
figurations analyzed have been either a laminate with a through-width or embedded 
delamination (see fig. 1). In this paper the embedded delamination will be 
considered. 

For the through-width delamination, both 2-D finite element [ref. 5] and 
Rayleigh-Ritz beam [ref. 7] solutions have been obtained. To date, strain-energy 
release rate analysis of the embedded delamination has been limited to using plate 
analysis to calculate the average value of the total strain-energy release rate G 
along a delamination front. The total strain-energy release rate is simply the sum 
of the mode I, mode II and mode III components. Ideally, one would like to obtain 
the magnitudes of the individual components. But such is beyond the scope of plate 
analysis; expensive three-dimensional analysis would be required. The current study 
also used plate analysis, so only total G was calculated. However, in contrast to 
earlier work, the present study calculated distributions of G along the delamina- 
tion front. The distribution of G was calculated using a new virtual crack closure 
technique which is suitable for use with plate analysis. 

The objectives of this paper are to present the derivation of the new virtual 
crack closure technique, evaluate the accuracy of the technique, and finally to 
present the results of a limited parametric study of laminates with a postbuckled 
delamination. Although the new virtual crack closure technique is general, only 
homogeneous, isotropic laminates were analyzed in this paper. This was to eliminate 
the variation of flexural stiffness with orientation, which occurs even for quasi- 
isotropic laminates. This made it easier to identify the effect of geometrical 
parameters on G. 

In the following sections the new virtual crack closure technique will be de- 
rived first. Then the specimen configurations will be described. Next, the stress 
analyses will be discussed. Finally, the virtual crack closure technique will be 
evaluated and then used to calculate the distribution of G along the delamination 
front of several laminates with a postbuckled de lamination. 
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submatrices in laminate theory constitutive relations 

length of delamination in x- and y-directions , respectively, mm 

flexural stiffness coefficient = EhV (1 2(1 -v 2 ) ) , Nm 

strain in x-direction in base laminate 

e Q for buckling of sublaminate 

2 

Young's modulus, N/m 

2 

total strain-energy release rate, J/m 

2 

average strain-energy release rate, J/m 

2 

strain-energy release rate due to moment, J/m 

2 

strain-energy release rate due to in-plane stress resultants, J/m 
curvature 

moment resultants, Nm/m 
in-plane stress resultants, N/m 
point load, N 

displacements in x, y, and z directions, m 

transverse displacement at the center of the square laminated 
plate, m 

rectangular Cartesian coordinates 
strain 

Poisson's ratio 


ANALYSIS 

The following subsections describe the new virtual crack closure technique 
( VCCT ) , the configurations analyzed, and the stress analyses used in conjunction with 
the VCCT. 


New Virtual Crack Closure Technique 

The basic tenet of any crack closure technique is that strain-energy release 
rate is equal to the work per unit area required to close a crack over a small dis- 
tance. For a 3-D configuration, such as that in figure 2, the strain-energy release 
rate is calculated using an equation of the following form [ref. 9]. 
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where Aa is the crack closure length, u, v, and w are the displacements in 
the x, y, and z directions, and 0 z* an< * a yz are ^ le tra ctions required 

on the crack faces to close the crack over t?ie distance Aa. The coordinate system 
is shown in figure 2. 

The various approximations made in plate theory make equation (1 ) not directly 
applicable for calculating G. For example, plate theory does not provide the magni- 
tude of the terms in equation (1). However, the concept of crack closure is still 
valid for plate analysis. The problem is to determine what is meant by closing a 
crack in plate analysis. In the following explanation of the new crack closure tech- 
nique, Kirchhof f -Love plate theory is assumed. 

Figure 3 shows the crack-front region of a plate with a single delamination. As 
shown in the figure, the intact plate (region A) is split at the crack-front into two 
sublaminates. Although the plate regions are shown with some thickness, the only 
gradient through the thickness is the linear variation of u and v due to the 
rotations -w, x anc * ” w 'y* Because of this simplified response in the thickness 
direction, there is an abrupt change in the response of the plate from the uncracked 
region to the cracked region. For example, the curvature k^ of region A at the 
crack front is different from k x of region B at the crack-front. Closure of the 
crack-front over an infinitesimal distance Aa implies the following: 


Over the distance Aa, the two sublaminates (regions B and C) behave as a 
single intact plate after closure. That is, u,v,e , e , and e vary 

X Y ^J: 

linearly and w is constant through the entire thickness of the combined 
plate (B + C) over the distance Aa. Since Aa is infinitesimal, the 
strain distribution through the thickness of the sublaminates after 
closure is the same as in the intact plate at the crack front. 


Work is required to impose the closure, since midplane strains and curvatures must be 
changed. This required work per unit area of crack closure is equal to the strain- 
energy release rate. 

The first step in the derivation of a mathematical formula for strain-energy 
release rate is to express the strain distributions in regions A, B, and C. These 
are given in equations (2). 
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To make the strain distributions through the thickness of regions B and C the 
same as for region A requires the following strain increments: 
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(3) 


Note that equations (3) express the strains in terms of each of the sublaminate 
coordinate systems. Equations (3) show that in general there is a change in both the 
midplane strain and the curvature. For each sublaminate equations (3) can be 
expressed as 
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where z* is either z - z B or z - z c , depending on which region is being 
considered. 


In equations (4) the twist curvature term Ak xy is zero. This is a consequence 

of the slope compatibility requirement at the crack* front, as explained next. 

Consider the crack front to be along a line x = constant. Along the crack front 

d w 9 w 

^ is continuous because of slope compatibility requirements. Since ^ is 

continuous regardless of the value of y , the change in in the y-direction. 
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To impose the required mid-plane strain and curvature changes Ae^ and Ak^ 
requires in-plane stress and moment resultants for each sublaminate. 

The magnitudes of the in-plane stress resultants and moments are determined 
using the sublaminate constitutive matrix, that is 
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The strain-energy release rate is simply the work required to change the mid- 
plane strains and curvatures in the cracked part of the plate at the crack front so 

that they are equal to the midplane strains and curvatures in the uncracked part of 

the plate at the crack front. This work is given by equation (7). 


0 - i ( 4 »x le xO * a »y 4e y° + 

q 3 B B \ 

+ AM Ak + AM Ak + similar terms for plate C ) 
y y xy xy / 


( 6 ) 


For a particular problem, many of the terms in equation (6) may be zero. If there is 
more than one delamination, the procedure is still the same, one simply has more sub- 
laminates to consider. Equation (6) will be used later in this paper to calculate 
the distribution of G for several configurations. 
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Cases Analyzed 


Two cases were examined: a transversely loaded square laminate and a post- 

buckled laminate in compression (fig. 4a). Figure 4b shows the plan view of the 
laminate and the delaminated region. The deformed profile is shown in figure 4c. 

The transversely loaded laminate was used to illustrate and evaluate the VCCT. 

It was selected because an exact solution is available for the displacements and 
boundary stress resultants [ref. 10]. This case consists of two isotropic plates 
bonded together except for a square region. The lower and upper plates will be 
referred to as the base laminate and the sublaminate, respectively. The base lami- 
nate is assumed to be rigid. The sublaminate has a Young's modulus of 207 GPa and a 
Poisson’s ratio of 0.3. The loading consists of a single transverse point load 
applied at the center of the sublaminate. The displacements are assumed to be small, 
so that linear analysis is adequate. 

The second case was selected to represent a compressively loaded laminate with 
an embedded delamination (see fig. 4a). The configuration is very similar to the 
transversely loaded laminate. The loading is in-plane (instead of transverse) and, 
although the lower plate is much stiffer than the top plate, it is not rigid. The 
sublaminate is isotropic with a Young’s modulus of 53.3 GPa and a Poisson’s ratio 
of 0.31. The thickness of the sublaminate is 0.51 mm. Again, the lower and top 
plates will be referred to as the base laminate and sublaminate, respectively. 

The deformation in the base laminate is assumed to be unaffected by the response 
of the sublaminate. Hence, the base laminate is in a state of uniaxial stress. The 
displacements of the base laminate are given by 


u = xe x 
v = -y v xy e x 


The slope at the sublaminate boundary is zero. 


Closed Form Analyses 

Two closed form stress analyses were used to analyse the transversely loaded 
square plate: a Rayleigh-Ritz (R-R) analysis [ref. 10] and an exact series solution 

[ref. 11]. The R-R analysis was used to evaluate the performance of the VCCT when an 
approximate stress analysis is used to obtain the terms in equation (7). The exact 
solution provided a reference for evaluating the accuracy of the VCCT when the terms 
in equation (7) are known exactly. The key equations from the R-R analysis in 
reference 10 are: 


\ \ [’ ' cos(ifli) . (. 

m=1 n=1 L N 


> n cos(^a) 


(7) 


6 



where 


a n = 0.12662 P 
a 22 = 0.00301 P 
a 33 = 0.0015 P 


2 a 2 -| a -0.00601 P 

a 1 3 = a 31 = 0.00278 P 


and the normalized load 
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Pa 

w 4 D 
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The maximum displacement w q under the load 


w Q = .005435 Pa^D^ 
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The key equations from the exact solution in reference 11 are 
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and an expression for the curvature k along x = a/2 
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-0.1025 P 
0.0263 P 
0.0042 P 
0.0015 P 


Only the first four terms were retained, since addition of the fifth term changed 
k x by only 0.32 percent. Because of symmetry, the distribution of k along 
y = a/2 is identical to that of k^ along x = a/2. Y 


Finite Element Analysis 

The STAGSOI [ref. 12] finite element program was used to analyse the con- 
figurations considered. Linear analysis was used for the transversely loaded 
laminate and geometrically nonlinear analysis was used for the postbuckled sub- 
laminate. The following subsections describe the finite element models and how the 
VCCT is used with finite elements. 
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Finite element models *- Figures 5 and 6 show typical finite element models used 
for the transversely loaded laminate and the postbuckled sublaminate, respectively. 
Boundary conditions are indicated in the figures. Because of symmetries about the 
x- and y-axis, only one fourth of the subplate is modeled for both configurations. 

The elements are non-conforming, four-node quadrilaterals with six degrees of freedom 
(three displacements and three rotations) at each node. The STAGSC-1 program refers 
to this element as type 410. The penalty function option was selected to reduce 
incompatibilities in the displacements at the element boundaries. 

The results presented herein were obtained with meshes graded in the x- and 
y-directions . This was found to give better results than a uniform mesh with the 
same number of elements. The grading of the mesh resulted in a fine mesh at one 
corner (eg. the upper right hand corner of the mesh in fig. 5). This fine mesh at 
the corner is just an artifact of the grading, not any' intention to model the corner 
more accurately than elsewhere. Several different mesh refinements were examined to 
assure that the trends were correct. 

Virtual crack closure technique for finite elements .- Theoretically, use of VCCT 
with finite elements is straight forward. Simply use the mid-plane strains and 
curvatures in regions A, B, and C (see fig. 3) at the crack front in equations (3), 
(4), (6), and (7). Unfortunately, many finite elements programs (including STAGSC-1) 
do not output the required quantities right at the crack front. For example, results 
may be provided at the quadrature points. In this case, one might extrapolate to the 
crack front. In general, some engineering judgement will be required to use the 
information provided by a finite element program to estimate the information required 
to calculate G. 

In this study, reactions at the boundary of the finite element model were used 
to calculate the required information. This was possible because only the sub- 
laminate was modeled; hence, the crack front was at the model boundary. The proce- 
dure used to extract the required information from the nodal reactions is described 
next. 


The boundary nodal reactions are the equivalent forces and concentrated moments 
which are required to impose the specified boundary conditions. For example the 
postbuckled sublaminate has specified nonzero in-plane displacements and zero slope 
along the edges x = ±a/2 and y = ±b/2. To impose these conditions requires a dis- 
tribution of mid-plane stresses and moments along these edges. The concentrated 
nodal moments and forces are statically equivalent to these distributed stresses and 
moments. A step variation of the moments was assumed along the edge (see fig. 7). 

The concentrated moment at a node was assumed to be equal to the integral of the 
distributed moment from the middle of the element on the left of the node to the 
middle of the element on the right. Nodes at the corners of the models had contribu- 
tions from only one side. This simple assumption means that the distributed moment 
is equal to the concentrated moment divided by the length. The mid-plane stress 
resultants were calculated using this same procedure. 
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RESULTS AND DISCUSSION 


First, details of the application and accuracy of the VCCT will be discussed for 
the transversely loaded laminate. The transversely loaded laminate is used because 
there is an exact, closed form solution for the plate deformations. Hence, the exact 
value of the average G along the delamination front can be calculated and used as a 
reference. Of course, it would be preferable to have a reference solution for the 
distribution of G, but one is not available. 

• 

Next, the effect of several parameters on G for postbuckled delaminations will 
% be discussed. 


Application of the VCCT 

Because of symmetry only one-eight of the perimeter of the transversely loaded 
plate need to be considered when applying the VCCT. In particular, the crack front 
along x = a/2 for 0 < y < a/2 will be analyzed. The plate deflections are as- 
sumed to be small. Hence, the mid-plane strains are zero throughout; therefore, 

Ae x o ~ ^ e y0 = ^ e xyO = Along x = a/2, the curvature k^ is zero for both the 
intact and the debonded region at the crack tip? therefore, Ak^ = 0. The curvature 
k x is zero throughout the lower plate, since it is assumed to be rigid. The top 
plate has zero k x in the bonded region and a non-zero k x in the debonded region. 
Therefore, the strain-energy release rate, using equation (7), is 


G - I (Ak M ) (11) 

2 xx 

where Ak * - (k of the unbonded top plate at the crack tip) 

X X 



and 


M x - D 1 1 Ak x 


Figure 8 shows the distribution of normalized G along x = a/2 obtained using 
several analyses. In all cases the VCCT was used to calculate G. The R-R solution 
used equation (7) to determine k x . The exact solution used equation (10) to 
determine k x * All three analyses predict a strong gradient of G along x = a/2. 
Growth of the delamination should occur preferentially in the middle of the sides, 
thus transforming the square into a rounded shape. The finite element results agree 
well with the exact solution. The R-R analysis predicts the same trend as the exact 
solution but the magnitudes are different. The large differences between the R-R and 
exact analyses are due to the errors in calculating boundary curvature and moment 
with the R-R analysis. 
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Average values of G for the crack front were calculated using the VCCT and by 
differentiating the strain energy with respect to M a". Although not a complete 
check, comparison of the results should give some indication of the accuracy of the 
VCCT. The average G was calculated using the VCCT by integrating the results in 
figure 8 and dividing by the length of the crack front (eq. (12)). 
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The strain-energy derivative technique (SEDT) uses the following expression for 
average G. 
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The derivative is taken with respect to a/2 because the delamination is assumed to 
extend in a self-similar manner. For example, the boundaries at x = + a/2 moves 
to x - ± (a/2 + A (a/2) ) • 

The strain energy U can be expressed as 



(14) 


From equations (13) and (14) we obtain 


G - 
av ' 4a 9a 


(15) 


The average G obtained using the VCCT and the strain-energy derivative technique 
(eq. (13)) are 

a av D 1l/ p2 


Exact 

0.002804 

(VCCT) 


0.002800 

(SEDT) 

Finite Elements 

0.002744 

(VCCT) 

R-R 

0.002229 

(VCCT) 


0.002717 

(SEDT) 
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There is essentially no difference between the average G obtained with the VCCT and 
equation (15) when the exact solution is used. (The small difference is a result of 
using only the first four terms in the series solution. ) The VCCT itself is appar- 
ently accurate, but the accuracy of the calculated G depends on the stress analysis 
used. When calculating average G with an approximate analysis, (i.e. the R-R anal- 
ysis) the VCCT is considerably less accurate than the strain energy derivative 
technique. 

Generally one would want to know the distribution of G rather than just an 
average value. If an approximate analysis is being used, equation (15) gives a bet- 
ter estimate of average G than using equations (11) and (12). A reasonable com- 
promise is to use equation (15) to obtain the average G and use equation (11) to 
obtain the distribution of G. Then scale the G distribution such that the VCCT 
gives the same average G. The curve labeled "scaled R-R" in figure 8 was obtained 
by using this procedure. The agreement with the exact solution is good. 


Postbuckled Sublaminate 

The VCCT described in the previous sections was used to calculate G for the 
compressi vely loaded specimen shown in figure 4. The effect of de lamination dimen- 
sions and strain level on the magnitude and distribution of G is illustrated in 
figures 9-12. 

For the postbuckled sublaminate, the equation for calculating G (eq. (6)) 
reduces to 


along x = a/2 

G = 1/2 ( AN 8 Ae® + AM® Ak® ) 

\ X xO X X / 


along y - b/2 

G = 1/2 ( AN® Ae® + AM® Ak® 

\ y yo y y 


(16) 
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Ae _ = e - e 
xO x 0 


Ae ^ » e + ve^ 
yO y 0 
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Ak = k 
x x 



e = strain in x-direction in base laminate 
o 



Equation (21) shows that G consists of contributions from changes in both mid-plane 
strain and curvature. The contributions from changes in the mid-plane strain and 
curvature are referred to herein as and G M , respectively. 

Figure 9 shows the distribution of G^ and G M for 25.4 x 25.4 mm and 
25.4 x 50.8 mm delaminations (solid and dashed curves, respectively). The first 
dimension refers to the length in the x-direction. The applied strains for the 
25.4 x 25.4 mm and 25.4 x 50.8 mm delaminations were -0.0040 and -0.0043 respec- 
tively. The bifurcation buckling strains were -0.0037 and -0.0018, respectively. 

For the square delamination (solid curves), the magnitudes of G^ and G M are 
much larger along y = b/2 (fig. 9b) than along x = a/2 (fig. 9a). Hence, one would 
expect an initially square delamination to grow transverse to the load direction (in 
the y-direction) . Since G M and G N are much smaller near the corners of the de- 
lamination, the delamination would not grow in a self-similar fashion. Instead, it 
would become more like an ellipse than a rectangle. Note that G^ is much larger 
than G n for the delamination front along both x = a/2 and y = b/2. For the 
rectangular (25.4 x 50.8) delamination the maximum value of the sum of G^ and G N 
is about 300 J/m 2 for both x = a/2 and y = b/2. Hence, one might expect the 
delamination to grow equally in the x- and y-directions, initially maintaining an 
aspect ratio (b/a) of about 2 and becoming elliptical. However, the aspect ratio for 
equal growth in the x- and y-directions depends on the absolute size of the delamina- 
tion and the strain level. This is discussed next. 

Figure 10 shows G (G = G {1 + G N ) for the 25.4 x 50.8 mm and 50.8 x 101.6 mm 
delaminations with an applied strain of -0.0034. Along y « b/2 (fig. 10b), the 
larger delamination has a G which is more than twice as large as G along x = a/2 
(fig. 10a). Hence, growth should occur more rapidly along y = b/2, resulting in an 
aspect ratio of more than 2. The smaller delamination in figure 10 was also sub- 
jected to a strain of -0.0034. Note that figure 10 predicts the smaller delamination 
will grow in the loading direction (because of the larger G), which will reduce the 
aspect ratio to less than 2. 

The peak values of G from figures 9 and 10 are combined in figure 11. The 
square delamination is clearly expected to grow more easily transverse to the load 
direction. The behavior of the 25.4 x 50.8 mm delamination depends on the strain 
level. At the lower strain level, growth is expected to occur preferentially in the 
load direction, but at the higher strain level, the growth rates are expected to be 
about the same in both directions. Although the 50.8 x 101.6 mm delamination has 
the same aspect ratio as the smaller rectangular delamination, figure 11 shows that 
the larger rectangular delamination will grow preferentially transverse to the load 
direction. Note that e 0 / e cR is 7.6 for the larger delamination, which is much 
larger than for the load cases for the smaller delamination. Apparently, as the 
strain becomes larger relative to the buckling strain, there is a greater tendency 
for the delamination to grow transverse to the load direction instead of along the 
load direction. 

Figures 9 and 10 suggest that at least for some configurations, delamination 
growth would be unstable if it is governed by G. In each case, the peak G is 
larger for the larger delamination when the applied strain is constant (approxi- 
mately). Of course, a rectangular delamination does not grow to be a larger rect- 
angular delamination; it becomes elliptical. But the large difference in G between 
the smaller and larger delaminations in figures 9 and 10 probably outweighs the 
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differences due to modeling as a rectangle instead of as an ellipse. More configura- 
tions need to be analyzed to determine whether the delamination arrests when it grows 
very large. 

The relative magnitudes of G N and G M depend both on the size of the de lami- 
nation and the edge examined. Figure 9 shows that for the 25.4 * 50.8 delamination 
G n is a significant protion of G along x = a/2, but it is insignificant along 
y = b/2. For the 25.4 x 25.4 mm delamination G N is insignificant along y = b/2. 
For the 25.4 x 25.4 mm delamination G N is insignificant along both edges for a 
strain of -0.004. 

Figure 12a shows the effect of strain level on the magnitudes of G N and G^ 
along x = a/2 for the 25.4 x 25.4 delamination. Dashed and solid lines are used 
for G m and G N , respectively. The peak value of G N decreased as the strain in- 
creased. Conversely, G^ increased monotonically . One of the G^ curves is marked 
"switched mode". For this strain level, the plate has changed from the first buck- 
ling mode to a higher mode. This can result in negative displacements, which are 
physically impossible since the base laminate will not allow this. Further work 
(beyond the scope of this paper) is needed to determine an accurate and efficient 
method for analyzing such cases. Figure 12b shows the effect of strain level on 
g m along y = b/2. The magnitude of G^ is not shown because it was quite small 
for all cases. The magnitude of G M increases rapidly as the strain increases. 


CONCLUSIONS 

A virtual crack closure technique was developed for calculating the distribution 
of total strain-energy release rate around the boundary of a delamination in a plate. 
This is an improvement over earlier techniques, which only provide an average of the 
strain-energy release rate along the delamination front. This technique was verified 
by analyzing a configuration with a closed-form exact solution for the average 
strain-energy release rate. The technique was then used to calculate the distribu- 
tion of strain-energy release rate for several cases involving the postbuckling of a 
delaminated region. The following conclusions were reached. 

1 • The virtual crack closure technique developed herein appears to be an ac- 
curate technique for calculating the distribution of total strain-energy release 
rate. 


2. For square and rectangular delaminations, there is a large variation of 
strain-energy release rate along the delamination front. Hence, self-similar growth 
is not expected. 

3. Whether a delamination grows in the load direction or perpendicular to the 
load direction in a compressively loaded laminate depends on the current delamination 
aspect ratio, the strain level, and the absolute size of the delamination. 
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(a) Deformed laminate and two types of loading considered. 



(b) Plan view of laminate. 



(c) Cross-section of deformed laminate. 


Figure 4.- Laminate configuration and loading 
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Figure 9.- Distribution of G^ and G N for two postbuckled 
sublaminate along x = a/2 and y = b/2. 
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Figure 10.- Distribution of G for two different size 
delaminations with same aspect ratio (e = 0.0034). 
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Figure 11.- Peak values of G for several sublaminate and loading cases. 
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Figure 12.- Effect of strain level on and for a = b = 

(G m is not shown in fig. 12(b), because it is very small for all 
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